﻿//---------------------------PROGRAM "LjuPad_1"----------------------------------
//
//-------------------ASTRONOMICAL OBSERVATORY OF PADUA--------------------------
//-----------------------AUTHOR: CLAUDIO FANTON---------------------------------
//
//            Ljubljana - Padova collaboration (Slovenia - Italy) 
//                   See Fanton et al, PASJ 1997, 49, 159
//                   Contact point: calvani@pd.astro.it
// Version 2.0, may  1998
//   major improvements:
//   - the program now uses an adaptive grid on the image plane to get a 
//     better image coverage. This gives smoother line profiles; 
//   - instead of binning the image area to calculate the line profiles, we now
//     evaluate the solid angle subtended by the disk in the observers' sky.  
// 
// Version 1.0, february 1998

/*
  This program has been tested on a Pentium PC (using Borland C++ 5.0 for Windows 95)
  and on a Dec Alpha OSF 433 Workstation using CXX.


  --------------------------------OUTPUT FILE FORMAT----------------------------


 The program "LjuPad_1" computes images for an accretion disk around a Kerr black
 hole.

 Each output data file corresponds to a particular image of a geometrically thin
 accretion disk. It has three lines of comment followed by the data:
  -  the first line describes the disk image parameters whose values
	 are written in the second line (all values are expressed in geometrized
	 units G=c=1; the radial coordinate "R" and the specific angular momentum
	 of the black hole "a" are adimensional quantities in units of M, the black
	 hole mass );
  -  the third line describes each column of the  file:
	  - the emission radius R;
	  - the g factor of the photon (that is the ratio
		between the observed and the rest frequency of the line);
	  - the solid angle subtended by the disk in the observers' sky;
	  - the two constants of motion (L and Q) of the photon trajectory.

 The accompaning program "LjuPad_2" that computes line profiles is adapted to
 read this particular file format.


  ----------------------------INPUT PARAMETERS--------------------------------

 The program will ask for the following input parameters
 (geometrized units c=G=1 are used):
 -  the name of the output file
 - "a": the specific angular momentum of the black hole in units of
		the black hole mass M;
 - "Theta_obs": the disk inclination angle to the line of sight.
		It is expressed in degrees. Theta_obs=0 is for a pole on accretion
		disk;
 - "R_in": the inner radius of the accretion disk in units of M.
		It must be >= R_ms, the radius of the marginally stable orbit;
 - "R_out": the outer radius of the disk in units of M;
 - "g_step": the step in g that will be used for binning the line profile.
	A good compromise is g=0.01, not to have too long computing times.
  ----------------------------------------------------------------------------
 */

#include <iostream>
#include <math.h>
#include <stdlib.h>
#include <fstream>
#include <iomanip>
#include <string>
#include <sstream>

using namespace std;

#define max(a, b)  (((a) > (b)) ? (a) : (b))

double a;
double val1, val2;

//metric functions
double delta(double r)
{
	return r * r + a * a - 2 * r;
}

double amet(double r, double th)
{
	return pow(r * r + a * a, 2) - a * a * delta(r) * pow(sin(th), 2);
}

double omega(double r, double th)
{
	return 2 * a * r / amet(r, th);
}

double kep(double r)
{
	return 1. / (pow(sqrt(r), 3) + a);
}

double sig(double r, double th)
{
	return r * r + a * a * pow(cos(th), 2);
}

double pph(double r, double lph)
{
	return r * r + a * a - a * lph;
}

double ro(double r, double fi, double th)
{
	return r * cos(th) / sqrt(1. - pow(sin(th) * cos(fi), 2));
}

double kk(double r, double th)
{
	return r * r / amet(r, th) * (r - 1. + delta(r) / amet(r, th) * (pow(r, 3) - a * a)) -
		(kep(r) - omega(r, th)) * (-1.5 * amet(r, th) * pow(kep(r), 2) / pow(r, 1.5) +
			(2 * a + (kep(r) + omega(r, th)) * (pow(r, 3) - a * a) / pow(r, 2)));
}

double ki(double r, double th)
{
	return sqrt(delta(r) * sig(r, th) / amet(r, th) - pow(sin(th), 2) *
		amet(r, th) / sig(r, th) * pow(kep(r) - omega(r, th), 2));
}


//computational functions
double amplitud(double uu, double emmc)
{
	double  ca = 0.0003, sn, cn, dn, aa, b, c, d, emc,
		u, u2, epu, emu;
	int l, i, ii;
	double em[13], en[13];
	int bo;

	emc = emmc;
	u = uu;
	if (emc != 0.0)
	{
		bo = (emc < 0.0);
		if (bo)
		{
			d = 1.0 - emc;
			emc = -emc / d;
			d = sqrt(d);
			u = d * u;
		}
		aa = 1.0;
		dn = 1.0;
		for (i = 0; i < 13; ++i)
		{
			l = i;
			em[i] = aa;
			emc = sqrt(emc);
			en[i] = emc;
			c = 0.5 * (aa + emc);
			if (fabs(aa - emc) <= (ca * aa))
				goto l1;
			emc = aa * emc;
			aa = c;
		}
	l1:   u = c * u;
		sn = sin(u);
		cn = cos(u);
		if (sn != 0.0)
		{
			aa = cn / sn;
			c = aa * c;
			for (ii = l; ii >= 0; --ii)
			{
				b = em[ii];
				aa = c * aa;
				c = dn * c;
				dn = (en[ii] + aa) / (b + aa);
				aa = c / b;
			}
			aa = 1.0 / sqrt(c * c + 1.0);
			if (sn < 0.0)
				sn = -aa;
			else
				sn = aa;
			cn = c * sn;
		}
		if (bo)
		{
			aa = dn;
			dn = cn;
			cn = aa;
			sn = sn / d;
		}
	}
	else
	{
		epu = exp(u);
		emu = 1.0 / epu;
		cn = 2.0 / (epu + emu);
		dn = cn;
		if (fabs(u) < 0.3)
		{
			u2 = u * u;
			sn = cn * u * (1.0 + u2 / 6 * (1 + u2 / 20 * (1 + u2 / 42 * (1 + u2 / 72))));
		}
		else
			sn = (epu - emu) / (epu + emu);
	}
	return sn / cn;
}

double compint(double qqc)
{
	double aa, aaa, b, bb, e, f, g, pp, em, t, q, ca = 0.0003, MPIO2;
	double qc;
	MPIO2 = 0.5 * 3.141592653589793;
	aa = 1.0;
	bb = 1.0;
	pp = 1.0;
	if (qqc == 0.0)
		cout << "Pause in routine ellint" << endl;
	qc = fabs(qqc);
	aaa = aa;
	b = bb;
	t = pp;
	e = qc;
	em = 1.0;
	if (t > 0.0)
	{
		t = sqrt(t);
		b = b / t;
	}
	else
	{
		f = qc * qc;
		q = 1.0 - f;
		g = 1.0 - t;
		f = f - t;
		q = q * (b - aaa * t);
		t = sqrt(f / g);
		aaa = (aaa - b) / g;
		b = -q / (g * g * t) + aaa * t;
	}
	while (1)
	{
		f = aaa;
		aaa = aaa + b / t;
		g = e / t;
		b = b + f * g;
		b = b + b;
		t = g + t;
		g = em;
		em = qc + em;
		if (fabs(g - qc) <= (g * ca))
			goto l3;
		qc = sqrt(e);
		qc = qc + qc;
		e = qc * em;
	}
l3: return MPIO2 * (b + aaa * em) / (em * (em + t));
}

double ellint(double x, double qqc)
{
	double ca = 0.0003, cb = 1.e-9;
	double qc, a, b, c, d, p, z, eye, y, f, em, e, g;
	int l;
	double ret;
	double aa = 1.0, bb = 1.0, MPI = 3.141592653589793;

	if (x == 0.)
		ret = 0.0;
	else
		if (qqc != 0.0) {
			qc = qqc;
			a = aa;
			b = bb;
			c = x * x;
			d = 1.0 + c;
			p = sqrt((1.0 + qc * qc * c) / d);
			d = x / d;
			c = d / (2.0 * p);
			z = a - b;
			eye = a;
			a = 0.5 * (b + a);
			y = fabs(1.0 / x);
			f = 0.0;
			l = 0;
			em = 1.0;
			qc = fabs(qc);
		_1:		b = eye * qc + b;
			e = em * qc;
			g = e / p;
			d = f * g + d;
			f = c;
			eye = a;
			p = g + p;
			c = 0.5 * (d / p + c);
			g = em;
			em = qc + em;
			a = 0.5 * (b / em + a);
			y = -e / y + y;
			if (y == 0.0)
				y = sqrt(e) * cb;
			if (fabs(g - qc) > ca * g) {
				qc = 2.0 * sqrt(e);
				l = l + l;
				if (y < 0.0)
					l += 1;
				goto _1;
			}
			if (y < 0.0)
				l += 1;
			e = (atan(em / y) + MPI * l) * a / em;
			if (x < 0.0)
				e = -e;
			ret = e + c * z;
		}
		else
			cerr << "Errore in Ellint" << endl;
	return ret;
}

void param_eq(double u1[4], double robs, double a, double qqph,
	double& coeff1, double& coeff2, double& coeff3,
	double& xrin)
{
	coeff2 = (2. * u1[2] - u1[0] - u1[3]) / (u1[0] - u1[3]);
	coeff3 = (u1[0] - u1[2]) * (u1[2] - u1[3]);
	coeff1 = 1. / sqrt(a * a * qqph * coeff3);
	xrin = coeff2 + 2. * (coeff3 + sqrt(coeff3 * (u1[0] - 1 / robs) * (1 / robs - u1[3]))) /
		((u1[0] - u1[3]) * (u1[2] - 1 / robs));
}

void param_mix(double u1[4], double u2[4], double robs, double a, double qqph,
	double& coeff1, double& coeff2, double& coeff4, double& coeff5,
	double& ttc, double& xrin)
{

	coeff4 = sqrt(pow(u1[1] - u1[0], 2) + pow(u2[1], 2));
	coeff5 = sqrt(pow(u1[1] - u1[3], 2) + pow(u2[1], 2));
	coeff2 = sqrt(coeff4 / coeff5 * (1 / robs - u1[3]) / (u1[0] - 1 / robs));
	coeff1 = 1. / sqrt(a * a * qqph * coeff4 * coeff5);
	ttc = sqrt((pow(coeff4 + coeff5, 2) - pow(u1[0] - u1[3], 2)) / (4. * coeff4 * coeff5));
	xrin = 2 * coeff2 / fabs(1 - coeff2 * coeff2);
}

void param_re(double u1[4], double robs, double a, double qqph,
	double& coeff1, double& ttc, double& xrin)
{

	coeff1 = 2. / sqrt(a * a * qqph * (u1[0] - u1[2]) * (u1[1] - u1[3]));
	ttc = sqrt((u1[0] - u1[3]) / (u1[0] - u1[2]) * (u1[1] - u1[2]) / (u1[1] - u1[3]));
	xrin = sqrt((u1[0] - u1[2]) / (u1[0] - u1[3]) * (1 / robs - u1[3]) / (u1[2] - 1 / robs));
}

int root_4(double a, double qqph, double lph, double* u1, double* u2)
{
	double ee, f, g, o, b, c, dddd, p, t, ddd, mod0, mod1;
	double y, x, x1, x2, x3, ovo, vov, ov, vo, absovo, absvov, MPIO2;
	double aa1, aa2, bb1, bb2, cc1, cc2, dd1, dd2, f0, f1;
	double s1[4], s2[4], r1[4], r2[4], mdl[4];
	double MPI = 3.141592653589793;
	int i, nu_re;

	MPIO2 = 0.5 * MPI;
	nu_re = 0;
	ee = -(lph * lph + qqph - a * a);
	f = 2 * (pow(lph - a, 2) + qqph);
	g = -a * a * qqph;
	o = 2 * ee;
	b = ee * ee - 4 * g;
	c = -f * f;
	if (f == 0)
	{
		dddd = ee * ee - 4 * g;
		if (dddd >= 0.0)
		{
			s1[0] = 0.5 * (-ee + sqrt(dddd));
			s1[1] = 0.5 * (-ee - sqrt(dddd));
			s2[0] = 0.0;
			s2[1] = 0.0;
		}
		else
		{
			s1[0] = -ee / 2;
			s1[1] = s1[0];
			s2[0] = 0.5 * sqrt(fabs(dddd));
			s2[1] = -s2[0];
		}
		mod0 = sqrt(pow(s1[0], 2) + pow(s2[0], 2));
		mod1 = sqrt(pow(s1[1], 2) + pow(s2[1], 2));
		if (s1[0] != 0.0 && s2[0] != 0.0)
			f0 = atan(s2[0] / s1[0]);
		if (s1[1] != 0.0 && s2[1] != 0.0)
			f1 = atan(s2[1] / s1[1]);
		if (s2[0] == 0.0 && s1[0] > 0.0)
			f0 = 0.0;
		if (s2[0] == 0.0 && s1[0] < 0.0)
			f0 = MPI;
		if (s1[0] == 0.0 && s2[0] > 0.0)
			f0 = MPIO2;
		if (s1[0] == 0.0 && s2[0] < 0.0)
			f0 = 3 * MPIO2;
		if (s2[1] == 0.0 && s1[1] > 0.0)
			f1 = 0.0;
		if (s2[1] == 0.0 && s1[1] < 0.0)
			f1 = MPI;
		if (s1[1] == 0.0 && s2[1] > 0.0)
			f1 = MPIO2;
		if (s1[1] == 0.0 && s2[1] < 0.0)
			f1 = 3 * MPIO2;

		r1[0] = sqrt(mod0) * cos(f0 / 2);
		r2[0] = sqrt(mod0) * sin(f0 / 2);
		mdl[0] = pow(r1[0], 2) + pow(r2[0], 2);
		r1[2] = -r1[0];
		r2[2] = -r2[0];
		r1[1] = sqrt(mod1) * cos(f1 / 2);
		r2[1] = sqrt(mod1) * sin(f1 / 2);
		mdl[1] = pow(r1[1], 2) + pow(r2[1], 2);
		r1[3] = -r1[1];
		r2[3] = -r2[1];
		u1[0] = r1[0] / mdl[0];
		u2[0] = -r2[0] / mdl[0];
		u1[2] = r1[2] / mdl[0];
		u2[2] = -r2[2] / mdl[0];
		u1[1] = r1[1] / mdl[1];
		u2[1] = -r2[1] / mdl[1];
		u1[3] = r1[3] / mdl[1];
		u2[3] = -r2[3] / mdl[1];
	}
	else
	{
		p = 1. / 3. * (3 * b - o * o);
		t = 1. / 27. * (2 * pow(o, 3) - 9 * o * b + 27 * c);
		ddd = pow(t / 2., 2) + pow(p / 3., 3);
		if (ddd > 0)
		{
			ovo = -t / 2. + sqrt(ddd);
			vov = -t / 2. - sqrt(ddd);
			y = 0.0;
			if (ovo != 0.0 && vov != 0.0)
			{
				absovo = pow(fabs(ovo), 1. / 3.);
				ov = ovo / fabs(ovo) * absovo;
				absvov = pow(fabs(vov), 1. / 3.);
				vo = vov / fabs(vov) * absvov;
				y = ov + vo;
			}
			if (vov == 0.0)
			{
				absovo = pow(fabs(ovo), 1. / 3.);
				ov = ovo / fabs(ovo) * absovo;
				y = ov;
			}
			if (ovo == 0.0)
			{
				absvov = pow(fabs(vov), 1. / 3.);
				vo = vov / fabs(vov) * absvov;
				y = vo;
			}
			x = y - o / 3.;
		}
		if (ddd == 0.0)
		{
			if (t != 0.0)
			{
				x1 = 2 * (-t) / fabs(-t) * pow(fabs(-t / 2), 1. / 3.) - o / 3;
				x2 = -(-t) / fabs(-t) * pow(fabs(-t / 2), 1. / 3.) - o / 3;
				x = max(x1, x2);

			}
			else
				x = -o / 3.;
		}
		if (ddd < 0.0)
		{
			x1 = 2 * sqrt(fabs(p) / 3) * cos(1. / 3. * acos((-t / 2) / sqrt(pow(fabs(p) / 3., 3))))
				- o / 3;
			x2 = 2 * sqrt(fabs(p) / 3) * cos(1. / 3. * acos((-t / 2) / sqrt(pow(fabs(p) / 3., 3)))
				+ 2 * MPI) - o / 3;
			x3 = 2 * sqrt(fabs(p) / 3) * cos(1. / 3. * acos((-t / 2) / sqrt(pow(fabs(p) / 3., 3)))
				+ 4 * MPI) - o / 3;
			x = max(x1, max(x2, x3));
		}

		//coeffients
		if (x >= 0.0)
		{
			cc1 = sqrt(x);
			cc2 = 0.0;
			aa1 = -cc1;
			aa2 = 0.0;
			bb1 = 0.5 * (ee + x + f / cc1);
			bb2 = 0.0;
			dd1 = 0.5 * (ee + x - f / cc1);
			dd2 = 0.0;
		}
		if (x < 0.0)
		{
			cc1 = 0.0;
			cc2 = sqrt(fabs(x));
			aa1 = 0.0;
			aa2 = -cc2;
			bb1 = 0.5 * (ee + x);
			bb2 = 0.5 * (f / cc2);
			dd1 = bb1;
			dd2 = -bb2;
		}
		s1[0] = pow(aa1, 2) - pow(aa2, 2) - 4 * bb1;
		s2[0] = -4 * bb2;
		s1[1] = pow(cc1, 2) - pow(cc2, 2) - 4 * dd1;
		s2[1] = -4 * dd2;
		mod0 = sqrt(pow(s1[0], 2) + pow(s2[0], 2));
		mod1 = sqrt(pow(s1[1], 2) + pow(s2[1], 2));
		if (s1[0] != 0.0 && s2[0] != 0.0)
			f0 = atan(s2[0] / s1[0]);
		if (s1[1] != 0.0 && s2[1] != 0.0)
			f1 = atan(s2[1] / s1[1]);
		if (s2[0] == 0.0 && s1[0] > 0.0)
			f0 = 0.0;
		if (s2[0] == 0.0 && s1[0] < 0.0)
			f0 = MPI;
		if (s1[0] == 0.0 && s2[0] > 0.0)
			f0 = MPIO2;
		if (s1[0] == 0.0 && s2[0] < 0.0)
			f0 = 3 * MPIO2;
		if (s2[1] == 0.0 && s1[1] > 0.0)
			f1 = 0.0;
		if (s2[1] == 0.0 && s1[1] < 0.0)
			f1 = MPI;
		if (s1[1] == 0.0 && s2[1] > 0.0)
			f1 = MPIO2;
		if (s1[1] == 0.0 && s2[1] < 0.0)
			f1 = 3 * MPIO2;
		s1[0] = sqrt(mod0) * cos(f0 / 2);
		s2[0] = sqrt(mod0) * sin(f0 / 2);
		s1[1] = sqrt(mod1) * cos(f1 / 2);
		s2[1] = sqrt(mod1) * sin(f1 / 2);
		r1[0] = 0.5 * (cc1 + s1[0]);
		r2[0] = 0.5 * (cc2 + s2[0]);
		mdl[0] = pow(r1[0], 2) + pow(r2[0], 2);
		r1[1] = 0.5 * (cc1 - s1[0]);
		r2[1] = 0.5 * (cc2 - s2[0]);
		mdl[1] = pow(r1[1], 2) + pow(r2[1], 2);
		r1[2] = 0.5 * (-cc1 + s1[1]);
		r2[2] = 0.5 * (-cc2 + s2[1]);
		mdl[2] = pow(r1[2], 2) + pow(r2[2], 2);
		r1[3] = 0.5 * (-cc1 - s1[1]);
		r2[3] = 0.5 * (-cc2 - s2[1]);
		mdl[3] = pow(r1[3], 2) + pow(r2[3], 2);
		u1[0] = r1[0] / mdl[0];
		u2[0] = -r2[0] / mdl[0];
		u1[1] = r1[1] / mdl[1];
		u2[1] = -r2[1] / mdl[1];
		u1[2] = r1[2] / mdl[2];
		u2[2] = -r2[2] / mdl[2];
		u1[3] = r1[3] / mdl[3];
		u2[3] = -r2[3] / mdl[3];
	}
	for (i = 0; i <= 3; i++)
	{
		if (u2[i] == 0.0)
			++nu_re;
	}
	return nu_re;
}

void sort_mix(double* u1, double* u2)
{
	double umax[4], ure[4], uim[4];
	int i, j, k, s;


	for (i = 0; i < 4; i++)
	{
		umax[i] = 0.0; ure[i] = 0.0; uim[i] = 0.0;
		if (u2[i] == 0.0)
			umax[i] = u1[i];
		else
		{
			ure[i] = u1[i]; uim[i] = u2[i];
		}
		u1[i] = 0.0; u2[i] = 0;
	}
	for (i = 0; i < 4; i++)
	{
		for (j = 0; j < (4 - i); j++)
		{
			if (umax[i + j] >= umax[i])
			{
				u1[i] = umax[i + j];
				umax[i + j] = umax[i];
				umax[i] = u1[i];
			}
		}
	}
	s = 0;
	for (k = 0; k < 4; k++)
	{
		if (umax[k] != 0)
		{
			u1[k] = umax[k]; u2[k] = 0;
		}
		else
		{
			for (i = s; i < 4; i++)
			{
				if (uim[i] != 0)
				{
					u1[k] = ure[i]; u2[k] = uim[i];
					s = i + 1;
					break;
				}
			}
		}
	}
}

void sort_re(double* u1, int n)
{
	double umax;
	int i, j;

	for (i = 0; i < n; i++)
	{
		for (j = 0; j < (n - i); j++)
		{
			if (u1[i + j] >= u1[i])
			{
				umax = u1[i + j];
				u1[i + j] = u1[i];
				u1[i] = umax;
			}
		}
	}
}

double test_eq(double rr1, double rr2, double coeff1, double coeff2,
	double xrin, double coeff3, double valpol, double* u1)
{
	double xco1, xco2, val1, val2;
	double r;
	xco1 = coeff2 + 2. * (coeff3 + sqrt(coeff3 * (u1[0] - rr1) * (rr1 - u1[3]))) /
		((u1[0] - u1[3]) * (u1[2] - rr1));
	xco2 = coeff2 + 2. * (coeff3 + sqrt(coeff3 * (u1[0] - rr2) * (rr2 - u1[3]))) /
		((u1[0] - u1[3]) * (u1[2] - rr2));
	val1 = coeff1 * (log(fabs(xco1)) - log(fabs(xrin)));
	val2 = coeff1 * (log(fabs(xco2)) - log(fabs(xrin)));
	val1 = fabs(val1 - valpol);
	val2 = fabs(val2 - valpol);
	if (val1 < 1.e-4)
		r = 1 / rr1;
	if (val2 < 1.e-4)
		r = 1 / rr2;
	return r;
}

double test_mix(double rr1, double rr2, double coeff1, double coeff4,
	double coeff5, double ttc, double xrin, double valpol,
	double* u1, int sgn_den)
{
	double co1, co2, co12, co22, xco1, xco2;
	double r = 0;


	val1 = 10;
	val2 = 10;
	if ((co12 = coeff4 / coeff5 * (rr1 - u1[3]) / (u1[0] - rr1)) > 0.0)
	{
		co1 = sqrt(co12);
		xco1 = 2. * co1 / fabs(1 - co1 * co1);
		if (sgn_den == 0)
		{

			val1 = coeff1 * (ellint(xco1, ttc) - ellint(xrin, ttc));

		}
		else
		{

			val1 = coeff1 * (ellint(xco1, ttc) + ellint(xrin, ttc) -
				2 * compint(ttc));

		}
	}



	if ((co22 = coeff4 / coeff5 * (rr2 - u1[3]) / (u1[0] - rr2)) > 0.0)
	{
		co2 = sqrt(co22);
		xco2 = 2. * co2 / fabs(1 - co2 * co2);
		if (sgn_den == 0)
		{

			val2 = coeff1 * (2 * compint(ttc) - ellint(xco2, ttc) -
				ellint(xrin, ttc));

		}
		else
		{

			val2 = coeff1 * (-ellint(xco2, ttc) + ellint(xrin, ttc));

		}
	}
	val1 = fabs(fabs(val1) - fabs(valpol));
	val2 = fabs(fabs(val2) - fabs(valpol));

	if (val1 < 1.e-7)
		r = 1. / rr1;
	if (val2 < 1.e-7)
		r = 1. / rr2;
	return r;
}



int main(int argc, char* argv[])
{
	double r, th, thobs, robs, rem;
	double z1, z2, rms, rin, rout;
	double rstep, gstep, fi, fistep;
	double x, y;
	double lph, qqph, g = 10;
	double mu2p, mu2m, qqc, xin, valpol;
	double xrin, ttc, ttc2, valrad, valrad1;
	double coeff1, coeff2, coeff3, coeff4, coeff5;
	double co1, co2, co3;
	double rr1 = 0.1, rr2 = 0.1, MPI = 3.141592653589793;
	double MPIO2;
	double u1[4], u2[4];
	int nu_re, sgnpolph, sgn_den = 0;
	int factor = 1;
	double f, sin_a, sin_b;
	double rph, rh, J;
	string nf;
	ofstream os;
	MPIO2 = 0.5 * MPI;

	//需要输入的参数顺序:黑洞自旋 a，黑洞倾角thobs，冷气体环内外半径rin，rout，红移因子gstep,输出文件名nf
	
	a = stod(argv[1]);
	thobs = stod(argv[2]);
	rin = stod(argv[3]);
	rout = stod(argv[4]);
	gstep = stod(argv[5]);

	//判断输入参数是否异常

	if (!(a >= 0.01 && a <= 0.9981))
	{
		cout << "Warning!!! Input a (rotational parameter a>=0.01 a<=0.9981)?: ";
		exit(1);
	}
	if (!(thobs >= 0.1 && thobs <= 89.9))
	{
		cout << "Warning!!! Input Theta_obs (Theta_obs>=0.1 Theta_obs<=89.9)?: ";
		exit(1);
	}
	z1 = 1 + pow(1 - a * a, 1.0 / 3.0) * (pow(1 + a, 1.0 / 3.0) + pow(1 - a, 1.0 / 3.0));
	z2 = sqrt(3 * a * a + z1 * z1);
	rms = 3 + z2 - sqrt((3 - z1) * (3 + z1 + 2 * z2)) + 0.001;
	if (rin < rms)
	{
		cout << "Warning!!! Input R_in (a=" << a << ", R_in>=R_ms=" << rms << ")?: ";
		exit(1);
	}
	if (rout < rin)
	{
		cout << "Warning!!! Input R_out (R_out>R_in=" << rin << ")?: ";
		exit(1);
	}
	if (gstep < 0.0)
	{
		cout << "Warning!!! Input g_step (g_step>0.00)?: ";
		exit(1);
	}
	os.open(argv[6], ios::out);
	if (!os)
	{
		cerr << "Couldn't open file: " << argv[6];
		exit(1);
	};

	thobs = thobs / 180 * MPI;    //angle in radiants

	rh = 1 + sqrt(1 - a * a);
	rph = 2 * (1 + cos(2. / 3. * acos(-a)));



	int rint = max(rin - 2., .2);
	rstep = gstep;
	fistep = cos(thobs) * gstep;
	for (rem = rint; rem <= ((rout + 2.) / cos(thobs)); rem += rstep)
	{
		if (rem >= rms)
			rstep = 2.5 * gstep * (1 + (rem - rin) / 5.);   // rstep is here!
		for (fi = 1.e-4; fi <= 2 * MPI + 1.e-4; fi += fistep)
		{

			th = thobs;
			fistep = 2.5 * cos(th) * gstep * sqrt(1. - pow(sin(th) * cos(fi), 2));  // fistep is here

			r = robs = 5000;
			x = ro(rem, fi, th) * cos(fi);
			y = ro(rem, fi, th) * sin(fi);

			if ((pow(x / rin, 2) + pow(y / (cos(th) * rin), 2)) >= 0.99 &&
				(x * x + y * y) <= pow(rout + 2., 2))
			{
				sin_b = (x * sig(r, th) / sqrt(amet(r, th)) / sqrt(amet(r, th) / delta(r) *
					pow(1. + omega(r, th) * x * sin(th), 2) + y * y + pow(x + a * sin(th), 2)
					- pow(r * r + a * a + a * x * sin(th), 2) / delta(r)));
				sin_a = -x * sig(r, th) * sqrt(delta(r)) / (amet(r, th) * sin_b *
					(1. + omega(r, th) * x * sin(th)));
				f = sig(r, th) * sqrt(delta(r)) / (amet(r, th) * sin(th)) + omega(r, th) *
					sin_a * sin_b;
				lph = sin_a * sin_b / f;
				qqph = 1. / delta(r) * (pow(pph(r, lph), 2) - pow(sig(r, th), 2) * delta(r) /
					(amet(r, th) * pow(sin(th), 2)) * ((1. - sin_a * sin_a) / (f * f))) - pow(lph - a, 2);

				if (y <= 0.)
					sgnpolph = 1;
				else
					sgnpolph = -1;

				// elliptic integral computation
					  //polar motion
				if (qqph > 0)
				{
					mu2p = 0.5 / (a * a) * (sqrt(pow(lph * lph + qqph - a * a, 2) + 4. * a * a * qqph) -
						(lph * lph + qqph - a * a));
					if (fabs(thobs) < 0.3)
						mu2p += 1.e-8;
					mu2m = 0.5 / (a * a) * (sqrt(pow(lph * lph + qqph - a * a, 2) + 4. * a * a * qqph) +
						(lph * lph + qqph - a * a));
					qqc = sqrt(mu2m / (mu2p + mu2m));
					if (pow(cos(th), 2) < mu2p)
					{
						xin = sqrt(mu2p / pow(cos(th), 2) - 1.);
						if (sgnpolph > 0)
							valpol = 1. / (a * sqrt(mu2p + mu2m)) * (factor * compint(qqc) - ellint(xin, qqc));
						else
							valpol = 1. / (a * sqrt(mu2p + mu2m)) * (factor * compint(qqc) + ellint(xin, qqc));

						//radial motion
						nu_re = root_4(a, qqph, lph, u1, u2);
						if (nu_re == 4)
						{
							sort_re(u1, 4);
							if (u1[1] != u1[2])
							{
								param_re(u1, robs, a, qqph, coeff1, ttc, xrin);
								ttc2 = ttc * ttc;
								valrad = valpol / coeff1 + ellint(xrin, ttc);
								valrad = amplitud(valrad, ttc2);
								r = (u1[0] - u1[2] + pow(valrad, 2) * (u1[0] - u1[3])) / ((u1[0] - u1[2]) * u1[3] +
									pow(valrad, 2) * (u1[0] - u1[3]) * u1[2]);
							}
							if (u1[1] == u1[2] && y < 0.0)
							{
								param_eq(u1, robs, a, qqph, coeff1, coeff2, coeff3, xrin);
								valrad = valpol / coeff1 + log(fabs(xrin));
								valrad = exp(fabs(valrad));
								co1 = pow((valrad - coeff2) * (u1[0] - u1[3]), 2) + 4. * coeff3;
								co2 = u1[2] * pow((valrad - coeff2) * (u1[0] - u1[3]), 2) + 2. * coeff3 *
									(u1[0] + u1[3] - (valrad - coeff2) * (u1[0] - u1[3]));
								co3 = pow(u1[2] * (valrad - coeff2) * (u1[0] - u1[3]), 2) + 4. * coeff3 *
									(coeff3 + u1[0] * u1[3] - u1[2] * (valrad - coeff2) * (u1[0] - u1[3]));
								rr1 = (co2 + sqrt(co2 * co2 - co1 * co3)) / co1;
								rr2 = (co2 - sqrt(co2 * co2 - co1 * co3)) / co1;
								r = test_eq(rr1, rr2, coeff1, coeff2, xrin, coeff3, valpol, u1);
							}
						} // end if 4 real roots
						else
						{
							sort_mix(u1, u2);
							param_mix(u1, u2, robs, a, qqph, coeff1, coeff2, coeff4, coeff5,
								ttc, xrin);
							ttc2 = ttc * ttc;
							if (1. - coeff2 * coeff2 > 0.0)
							{
								valrad = valpol / coeff1 + ellint(xrin, ttc);
								valrad1 = 2 * compint(ttc) - ellint(xrin, ttc) - valpol / coeff1;
							}
							else
							{
								sgn_den = 1;
								valrad = 2. * compint(ttc) - ellint(xrin, ttc) + valpol / coeff1;
								valrad1 = ellint(xrin, ttc) - valpol / coeff1;
							}
							valrad = amplitud(valrad, ttc2);
							valrad1 = amplitud(valrad1, ttc2);
							if (valrad > 0.0)
							{
								co1 = (-1. + sqrt(1 + pow(valrad, 2))) / valrad;
								rr1 = (coeff4 / coeff5 * u1[3] + co1 * co1 * u1[0]) /
									(coeff4 / coeff5 + co1 * co1);
							}
							if (valrad1 > 0.0)
							{
								co2 = (1. + sqrt(1 + pow(valrad1, 2))) / valrad1;
								rr2 = (coeff4 / coeff5 * u1[3] + co2 * co2 * u1[0]) /
									(coeff4 / coeff5 + co2 * co2);
							}
							r = test_mix(rr1, rr2, coeff1, coeff4, coeff5, ttc, xrin, valpol,
								u1, sgn_den);
						} // end else 4 real roots
					}//end if cos(th)
				} // end if qqph
	// end elliptic integrals

				if (r >= rin && r <= rout)
				{
					th = MPIO2;
					g = 1. / (1. - kep(r) * lph) * ki(r, th);
					J = ro(rem, fi, thobs) / rem * rstep;
					J = J * fistep * (ro(rem, fi, thobs) + J / 2.);


					os << x << "  " << y << "  " << r << "  " << g << "  " << J << "  " << lph << "  " << qqph << endl;
					//               osx<<x<<'\n';
					//				osy<<y<<'\n';
					//				osx<<x<<"  "<<y<<"  "<<r<<endl;
					//				osy<<y<<endl;
					//				osr<<r<<endl;
				}
			}
		}
	}
	os.close();
	//   osx.close();
	//   osy.close();
	//   osr.close();
	return 0;
}
